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ABSTRACT 

This paper considers compressible turbulent flows at low turbulent Mach numbers. Con- 
trary to the general belief that such flows are almost incompressible, (i.e. the divergence 
of the velocity field remains small for all times), it is shown that even if the divergence of 
the initial velocity field is negligibly small, it can grow rapidly on a non-dimensional time 
scale which is the inverse of the fluctuating Mach number. An asymptotic theory which en- 
ables one to obtain a description of the flow in terms of its divergence-free and vorticity-free 
components has been developed to solve the initial- value problem. As a result, the various 
types of low Mach number turbulent regimes have been classified with respect to the initial 
conditions. Formulae are derived that accurately predict the level of compressibility after the 
initial transients have disappeared. These results are verified by extensive direct numerical 
simulations of isotropic turbulence. 
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No. NAS1-18605 while all of the authors were in residence at the Institute for Computer Applications in 
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I. Introduction 


Turbulence is the most common state of fluid motion. It is an all pervasive ubiq- 
uitous phenomenon present in, to cite a few instances, weather patterns, ocean cur- 
rents, high-temperature plasmas, astrophysical jets, and combustion. Despite the 
best theoretical and experimental attempts of more than a century, and more recent 
computational approaches, turbulence remains to be one of the great unsolved prob- 
lems of fundamental physics, and poses a grand challenge to any type of scientific 
investigation. 

There has of course been some progress in our understanding of turbulence in 
low-speed flows which are essentially incompressible (see Dwoyer, Hussaini and Voigt 
1984). There have also been a few attempts at predicting the effect of compressibility 
on some known results of incompressible turbulence such as the Kolmogorov spectrum 
(Zakharov and Sagdeev 1970, Kadomtsev and Petvishvili 1973, Moiseev, Sagdeev, Tur 
«*nd Yanovskii 1977, and L’vov and Mikhailov 1978). The majority of work in com- 
pressible turbulence focuses naturally on the linear regime, owing to the inherent 
difficulty of treating nonlinearities which now involve variable density as well. Moyal 
(1952) appears to be the first one to look at spectra of homogeneous isotropic turbu- 
lence in compressible fluids. He decomposed the velocity field in the spectral space 
into a longitudinal component (random noise) and a transverse (eddy turbulence) 
component. This is in fact equivalent to a Helmholtz decomposition (gradient of a 
scalar and curl of a vector) in physical space. The longitudinal component in physical 
space is variously known as the acoustic component or compression turbulence; the 
transverse component is also termed the incompressible component, the solenoidal 
component or shear turbulence. A broad conclusion of his work was that the inter- 
action between these two components is exclusively due to the nonlinear terms, and 
such interactions are the strongest at high levels of turbulence and at high values of 
Reynolds number. The properties of sound emitted in the far field by eddy turbulence 
were studied by Lighthill (1952, 1954, 1956). His formula for sound emission is found 
in laboratory experiments to remain valid far beyond the linear regime. Kovasz- 
nay (1957) proposed a decomposition of compressible turbulence into three modes - 
the vorticity mode, the entropy mode and the acoustic mode — and showed how to 
determine their levels and correlations from mass flow and stagnation temperature 
fluctuations measured by a hot-wire anemometer. Chu and Kovasznay (1958) have 
outlined a consistent successive approximation procedure in terms of an amplitude 
parameter, and have provided explicit formulae for second-order interactions among 
these three modes. They provide no explicit solutions since their main purpose was 
to provide a consistent framework to assess the nonlinear interactions in the experi- 
mental data. Tatsumi and Tokunaga (1974) and Tokunaga and Tatsumi (1975) study 
the interactions of weakly nonlinear disturbances such as compression waves, expan- 
sion waves and contact discontinuities using a reductive perturbation method due 
to Taniuti and Wei (1968). The key result is that the interaction between waves 
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of different families of characteristics leads to alterations in their amplitudes, phase 
velocities and propagation directions, whereas the interaction between waves of the 
same family of characteristics causes merger or coalescence. They further inferred that 
the statistical properties of two-dimensional shock turbulence are similar to those of 
one-dimensional shock turbulence which in turn are identical to those of Burgers tur- 
bulence. More recently, in a preliminary study of return to isotropy in a compressible 
flow, Lumley (1989) has provided some orders of magnitude estimates in terms of 
Mach number and Reynolds number. An excellent review of the published literature 
on compressible turbulence up to 1967 may be found in Monin and Yaglom (1967). 

The computational approach to turbulence is based on the Navier-Stokes equa- 
tions. The first attempt to solve numerically the equations of motion for compressible 
homogeneous turbulence is due to Feiereisen, Reynolds and Ferziger (1981). They as- 
sumed the divergence of the initial flow field and its time-derivative to be both zero. 
It was therefore not surprising that their results for isotropic compressible turbulence 
did not show any departure from the corresponding incompressible data. Recently, 
Passot and Pouquet (1986) have carried out numerical simulations of two-dimensional 
homogeneous compressible turbulent flows. They show that the behavior of the flow 
beyond an initial turbulent Mach number of 0.3 differs sharply from the lower Mach 
number cases which are characterized by the absence of shocks. 

In the present paper, we develop an asymptotic theory (with turbulent Mach 
number as a small parameter) to solve the compressible Navier-Stokes equations. 
The problem is set up as an initial- value problem to study the influence of initial 
conditions on the subsequent turbulence structure and its dynamical evolution. Ex- 
plicit relationships between the initial turbulent Mach number, pressure and velocity 
fluctuation levels are derived which are valid after the initial time transients have 
disappeared. Direct numerical simulations of two-dimensional isotropic compressible 
turbulence are performed to validate the results. 


II. Theoretical Analysis 


A . Problem Formulation 


The compressible Navier Stokes equations are, in non-dimensional variables, 

0 , 




dt 


+ V-(pun) 


dP 

— — (- u • VP + 7 -PV-u 
dt 


- 7 MA VP+ Pe V<r ’ 

7 V-kVT + 7(7 ~ 1)M ^, 


( 1 ) 
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where 


a = 2fi 


L i(Vu + Vu T )-^(V-u)/ 


( 2 ) 


is the viscous stress tensor (with bulk viscosity assumed zero), I is the identity tensor, 
and 

$ = -(Vu + Vu T ) : a (3) 

2 

is the viscous dissipation function. The equation of state is 

P = pT. (4) 


The density, velocity, temperature and pressure are respectively 

Pi 

PR 

u* 

Ur 

II 

Tr ’ 

P* 

R 3 PrTr 

where the dimensional reference values are denoted by the subscript r and R g is the 
universal gas constant. Dimensional variables have an asterisk superscript. Distance 
and time are scaled respectively with respect to Lr and tR = Lr/Ur. In the text, x 
refers to the cartesian position vector. 

The Prandtl number Pr = (prC p )/ kr, the Reynolds number Re — ( PrUrLr )/ fiR\ 
7 is the ratio of specific heats. Viscosity and conductivity are denoted by p and k 
respectively and are scaled with respect to hr and kr. 

The reference Mach number 

( 5 ) 

is related to the time-dependent turbulent Mach number 



P = 

u = 

T = 
P = 


M t 



u 


•2 


7 R 3 T* 


> 


Pt> 


( 6 ) 


The objectives of the following analysis is first to classify the various types of 
equilibrium turbulent regimes (distinguished by the presence or absence of shocks and 
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also by the fraction of the total kinetic energy solely due to compressibility effects), 
and second, to predict the range of initial conditions that leads to each one. The 
effect of viscosity is felt either on a viscous time scale (much greater than the acoustic 
time scale), or during the formation of shocks. In the latter case, viscosity serves to 
prevent the formation of singularities. Although there is a distinct possibility that 
shocks will form as a result of certain types of initial conditions, viscosity does not help 
initiate the processes (wave steepening) which eventually lead to shocks. Rather, the 
viscosity diffuses the sharp gradients to generate shocks of finite thickness. The above 
considerations lead us to neglect all viscous effects in the theoretical formulation of 
the full viscous, Navier-Stokes equations. This assumption will be verified a posteriori 
by the direct numerical simulations. 

In the following analysis, the turbulent Mach number is assumed to be very much 
less than unity, so that the sound velocity is much greater than the flow velocity. 
Under these circumstances, the acoustic time scale is much less than the convective 
time scale, which in turn is much less than the viscous time scale (if the Reynolds 
number is sufficiently large). As will be shown later, the different regimes emerge on 
the acoustic time scale, during which time any inconsistency in the acoustic compo- 
nent of the flow is washed away. In other words, time boundary layers only extend 
over a time period of 0(M R ). 

Dropping the viscous and heat conduction terms, Eqs. (1) become (in nonconser- 
vative form) 

dp 

— + u • Vp + pVu = 0, 

, du 1 

K _ + u .V u ) + — VP = °, (7) 

dP 

— + u • VP = — yPV-u. 

The initial conditions are arbitrary with the restriction that both Mr and the root 
mean square (rms) pressure, ||p(x, 0)||, be much less than unity. For the remainder of 
this analysis, we assume that for homogeneous turbulent flows, the maximum norm 
and the L 2 norm are of the same order of magnitude. 

B. Asymptotic Analysis 

The system of equations (7) is now investigated by an asymptotic analysis which 
assumes M R « 1. The pressure P is decomposed into a mean and a fluctuating 
component: 

p = 1+ P, (8) 

where 

I Ml « 1- (9) 
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For simplicity, we neglect any consideration of the entropy mode and assume the 
homentropic condition 

P = p\ (10) 

The density is expanded in powers of p according to 

P = l + 7 _1 P+0(||p|| 2 ). (11) 

Using the expansions (8) for pressure and (11) for density, Eqs. (7) become (to first 
order in ||p||) 


du 1 

aT + uVu + ^ | Vp = °- 


dp 

-^7 + u • Vp 4- 7 V -u = 0. 
ot 


( 12 ) 


At t = 0 the initial data is 


u(x,0) = u 0 (x), 
p(x,0) = p 0 (x). 


(13) 


It is convenient for the purpose of analysis to split the velocity vector according 
to 

u = u 7 + u c , (14) 

where the solenoidal component u 7 (x, t) satisfies the time-dependent incompressible 
Navier-Stokes equations 


du 1 

~dt 


+ u 1 ■ Vu 7 


-Vp 7 , 


Vu 7 = 0, 

with initial conditions 

u 7 (x,0) = u 7 (x), 

V-u 7 = 0. 

The incompressible pressure p 7 satisfies the Poisson equation 

V 2 p 7 = -Vu 7 : Vu 7 

from which it is inferred that 

p 7 = 0(u 72 ). 


(15) 


(16) 


(17) 

(18) 
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The initial solenoidal velocity is determined from the decomposition 

Uo(x) = u 7 (x) + u£(x). (19) 

If the flow is homogeneous, it is easy to show that uj(x) and Uq(x) are unique 
once Uo(x) is specified. When the flow is non-homogeneous, an arbitrary potential 
function can be added to the decomposition (19). The transverse velocity component, 
u 7 is divergence- free for all time (since it satisfies the time-dependent incompressible 
Navier-Stokes equations). Although u c is initially vorticity-free, it acquires a small 
degree of vorticity subsequently. Kreiss, Lorenz and Naughton (1990) have estimated 
this vorticity to be 0(Mj i ). 

Substitution of Eq. (14) into Eqs. (12) and manipulation of the resulting expres- 
sions using Eqs. (15) yield 

+ u 7 • Vu c + u c • Vu ; + u c • Vu c = -^y(Vp- 7M^Vp 7 ), 

( 20 ) 

757 + u 7 • Vp + u c • Vp + 7V u c = 0. 
ot 

The momentum equation in Eqs. (20) is simplified if the perturbation pressure is 
decomposed according to 

p = 7 M*p 7 + Sp c . (21) 

This particular decomposition removes p 7 from the evolution equation for u c . At 

t = 0, 

p 7 (x, 0) = p 7 0 (x) (22) 

is uniquely determined from the incompressible evolution equations, so that there 
is no ambiguity in the computation of p c (x,0). The uniqueness of the pressure 
decomposition is maintained for all time. The parameter 8 is defined to insure that 

||p c (x,0)|| = 1. (23) 

The total initial velocity u is also normalized to unity, i.e. ||u 0 (x)|| = 1. This is most 
easily accomplished by choosing 

Ur= |K(x)||. (24) 

Therefore, both u£ and u£ are 0(1). From Eq. (18), ||p£|| = 0(1). Incompressible 
variables vary on a slower time scale than the acoustic time scale, which is the time 
scale of interest, so that their order of magnitude will remain invariant in time. In 
contrast, we will demonstrate that certain combinations of initial conditions can lead 
to time boundary layers which can change the order of magnitudes of compressible 
variables on a time scale of O(Mr). 
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Substitution ofEqs. (14), (21) into the system (12) leads to the evolution equations 

= 0 , 


+ u 7 ■ Vu c + u c • Vu c + Au c -\ tttVp c 


dt 


dp c 


7 ma 

+ u 7 • Vp c + u c • Vp c + Bu c + yV-u c = G, 


dt * r 6 

for the compressible components of the flow variables, where 

Au c = u c • Vu 7 , 

-u c • Vp 7 , 


Bu c = - ,M *- 


G = 


— ( aT + " ' Vp 1- 


For reference, the orders of magnitude for A, B and G at t = 0 are 

•4 = 0 ( 1 ), 

B = 


G = 0( 




Since A, B, G only depend on incompressible quantities, their order of magnitude i 
constant over the time scale of interest. The initial data is 


is 


u c (x,0) = U 0 (x)-u 7 (x) = u£(x), 

p c (x,0) = Po(x)- ^p£ = Po(x). 

Note that is vorticity-free. The pressure -P(x) is now given by 

P(x, t) = 1 + 7 M£p 7 (x, t ) + 6p c ( x , t). (25) 

It must be noted that although pP and are 0(1) at t = 0, there is no guaranty 
that they will remain so. Should they exceed that order, they will be scaled down, 
otherwise terms might be neglected in the governing equations which should be kept. 


G. Regime Classification 


Freezing coefficients in Eqs. (25) makes it clear that there are two different time scales 
present, i.e. the convection scale of order 0(1), and the sound wave scale of order 
0(M r ). Since the initial compressible velocity field u£(x) is vorticity-free, we expect 
f° vary on the fast acoustic time scale. Therefore we introduce a new time variable 



( 26 ) 


7 



0, 


In terms of Eqs. (25) become 

^ + M R [(u I + u c )-Vu c + Au c } + ^—Vp c = 

( 27 ) 

+ Mr [( u 7 + u c ) • Vp c + Bu c ] + ^y^V-U C = Mr G. 

The system of equations (27) is hyperbolic and for 8 « Mr or 8 » Mr it is highly 
asymmetric. Therefore the solution will in general grow rapidly, thus invalidating 
the original scalings (see appendix). Only with the symmetrized version of Eqs. (27) 
can we be sure that the solution at t > 0 is bounded by the order of magnitude of 
the solution at t = 0. Therefore we symmetrize the system in such a way that the 
initial data are not “scaled up” because otherwise nonlinear terms which are initially 
negligible may not remain so. The change of variables which leads to symmetry is 
dependent on whether the coefficient of Vp c in Eq. (27) is less or greater than unity. 
In terms of 8 , we therefore distinguish two cases. 

8 < Mr') 


This regime is characterized by low levels of initial acoustic pressure fluctuation. 
We first introduce a new dependent variable 


which insures that 
Eqs. (27) become 


Mrt) 


Ill'll <||rf|| = 0(1). 


dp c ' 

dt 

with initial data 


d\i^ 

-Qp- + Mr[( u 1 + u c ) • Vu c + Au c ] + Vp c ' 

+ Mr [( u 7 + u c ) • Vp c H f--Bu c ] + V u c 

7 Mr 


0, 



(28) 

(29) 


(30) 


u c (x, 0) 

P C '(x,0) 


»?(*), 



By construction, u c (x, 0) and p c, (x, 0) are 0(1). If they should be much less than 
unity, the dependent variables are rescaled so that they become 0(1), in which 
case the coefficients of Eqs. (30) are replaced by smaller coefficients, and the ap- 
proximations below become more valid. (The notation / = O(M) is equivalent to 
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hm M o f ( M )/ M ^ 0 < oo. On the other hand, the O(M) symbol only states that 
the hmit is finite, possibly zero.) The system (30) is a well-behaved system of dif- 
ferential equations. However, one must check whether the convective terms remain 
negligible when t' = 0(1). It is easily shown that on the fast time scale, the lowest 
order evolution equations are given by the linear wave equations 

= o, 


= o, 

with the initial data 


C7U . 

dp c ' 


df 


+ V-u c 


u C ( x , 0) = 

P C '(x,0) = 

On the convective time scale, Kreiss et al. (1990) have shown that the wave equation 
still describes the acoustic component of the flow. 

From the wave equation, it is now easy to predict the order of magnitude of u c 
and p as a function of initial magnitudes. After a time t' = 0(1), both u c and p c ' 
will have the same order of magnitude, given by 

u c (x,0*J> c '(x,0 = 0(max(u£(x),p£'(x))), 

= 0(max(u£(x), —£—)). 

IMr" 

The results are summarized in table 1 for all different initial estimates of velocity. 


«?(*), 


u c o 

P C 0 

r « c » 

P u 

r oo 

0(1) 

O(M’r) 

0(1) 

0(1) 

0(1) 

0(max(M£, ^-)) 

0(max(^-,l)) 


Table 1: Order of magnitude of equilibrium levels of compressible pressure and 
acoustic velocity as a function of initial levels. 


c Wi ^ * he ctosen initial normalization for the velocity, u' must be 0(1), unless 
u — (1). The choice of the initial magnitude of the incompressible velocity compo- 

nent partly determines the initial variation of the ratio of compressible kinetic energy 
to total kinetic energy. 

£ > Mr 7 
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Now we introduce the new variable 


u c. = M*I u c 


which insures that 

iiug'ii < Hug'll = o(i). 

Using Eq. (31), Eqs. (27) become 

du c ' ... r 8 


dt 


— + (Mru 1 + -u c ') ■ Vu c ' + MrAvF' + Vp' 
7 7 


c _ 


= 0 , 


%- + {Mru 1 + -U c ') • Vp c + - B\ 1 C ' + V-u 1 
dt' K 7 7 


c f _ 


= A/k G 


with the initial data 


u c '(x,0) = ^-^Uo(x), 


p c (x,0) = Po(x) = 0(l). 


(31) 

(32) 


(33) 


(34) 


Once again, if the initial data in Eq. (34) should be much less than unity, a rescaling 
of velocity and pressure would simply decrease the coefficients of the quadratic terms. 
As long as u c '(x,f') = 0(1), Eqs. (33) reduce to 


du c> 

dt' 


+ 


-u c ' • Vu c ' + Vp c 
7 


0, 


dp c 

w 


+ 


-u c '- Vp c + -Bu c ' + V-u c ' = 
7 7 


0, 


(35) 


with initial conditions given by Eq. (34). When t' — 0(1), 


||u c '(x,i / )|| = ||p c (x,fOII, 

= 0(max(po,Uo')), 

= 0 ( 1 ). 


Equations (31) and (36) imply that 




(36) 


(37) 


When 8 = 0(1), the convective terms balance the time derivative terms, in which 
case wave steepening may occur on a t 1 — 0(1) time scale, and there is a propensity 
for shocks to form. Note that when shocks do occur, the turbulent Mach number, M t 
is no longer small, but has increased by a factor 1 /Mr from its initial value. 
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D. Results 


In the previous section, we established the conditions for two separate asymptotic 
regimes, depending on whether the ratio of 8 to ^M R was lesser or greater than unity. 
We now wish to deduce some physical consequences from the results of the previous 
section. To this end, we consider various levels of initial rms values for u c 0 and 6. 
The rms of u J is fixed at unity. Only when ||u£|| = 0(1) can ||u 7 || < 0(1). It is then 
easily shown that the rms level of can be increased to unity without affecting the 
order of magnitude estimates at t' = 1. 


Table 2 summarizes the theoretical findings from the previous section. To properly 
interpret the table, replace a cell value of, for example, 3 by 0(M R ). A subscript ()„, 
refers to the value of a variable after the initial transients have died away ( t ' = 0(1)), 
also referred to as an equilibrium state. When viscous effects are taken into account, 
the equilibrium state of all quantities varies over the viscouB time scale. They also 
fluctuate in time due to higher order convective effects that have been neglected to 
lowest order in the current expansion. For example, on the incompressible convective 
time scale, the terms Mru 1 • Vu c induce small fluctuations of the variables about 
the equilibrium state of the system. 

The rms levels of u£ and 6 are given in the first two columns of table 2. The fifth 
and sixth columns contain the new quantity 



which is the fraction of the total kinetic energy solely due to acoustic effects. Note 
that at t — 0, the denominator is unity by construction, so that 



and 

IKII = v/5S- 


(40) 


From Eq. (37), 



(41) 


and since ||u 7 || = 0(1), 8 = 0(1) always implies that x«> = 0(1). Consequently, 
the lower the initial level of compressibility at a fixed M R (i.e. the lower xo), the 
stronger the initial transient (which extends over the time period 0(M R )). As this 
imbalance intensifies, the propensity for shocks to form also increases. Table 2 re- 
flects these statements in the rows corresponding to 8 = 0(1). When Uq = 0(1) 
u' = O(Mr) j n > 1 . However, the value of n does not change the previous conclu- 
sions. Based on the results in Sarkar, Erlebacher, Hussaini and Kreiss (1989), it is 
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straightforward to show that 


X °° 2 1 -0.5xo(l -( 2 Y 

where „ 

-Mi 
5 M,j' 


(42) 


(43) 


u 0 

8 

U C ’oo 

P d oo 

Xo 

Xoo 

^0 

0 

0 

-1 

0 

0 

0 

2 

0 

1 

0 

0 

0 

0 

0 

0 

2 

0 

-1 

0 

0 

-2 

0 

3 

0 

-2 

0 

0 

-4 

l 

0 

-1 

0 

2 

0 

4 

l 

1 

0 

0 

2 

0 

2 

i 

2 

1 

0 

2 

2 

0 

l 

3 

1 

-1 

2 

2 

-2 

l 

4 

1 

-2 

2 

2 

-4 

2 

0 

-1 

0 

4 

0 

6 

2 

1 

0 

0 

4 

0 

4 

2 

2 

1 

0 

4 

2 

2 

2 

3 

2 

0 

4 

4 

0 

2 

4 

2 

-1 

4 

4 

-2 

2 

5 

2 

-2 

4 

4 

-4 


Table 2: Estimates of equilibrium levels of compressible pressure (p c ), (u c ), 

X as a function of a wide range of initial levels for p c and u c . Table entries 
should be interpreted as 0(M t ^ ble cntrv ). 

Two separate regimes (different from the mathematical regimes of the previous 
section) are apparent from table 2. For large values of 8 , ||p£j| = 0(1), while u c 
increases sharply. This sudden increase over the acoustic time scale is most important 
when 8 is large and Uq is small. The fact that this increase in the acoustic velocity is 
not contingent on 8 ~ 0(1) indicates- that an initial imbalance in compressible energy 
can occur without the need to appeal to wave steepening. When 8 is sufficiently 
small, u c remains at its initial level and large pressure waves appear in the flow. 
The compressible pressure imbalance becomes more severe with decreasing Uq and 
decreasing 8 . Equilibrium initial conditions are characterized by the juncture between 
the two aforementioned regimes. This corresponds to 

“0 = < 44 > 
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Using Eq. (39), Eq. (44) becomes 


= 6(1). (45) 

The function 

m = ^ ( 46 ) 

was introduced by Sarkar et al. (1989) and shown to play a fundamental role in the 
description of compressible turbulence. The last column in table 2 tabulates the value 
of Fq = F(0) for the various cases. It is now easy to see that a necessary condition 
for the appearance of shocks is 6 = 1 (see previous section). If the latter condition 
is not satisfied, the level of compressibility will remain small for all time, which is 
inconsistent with the appearance of an isotropic distribution of shocks. 

It is interesting to establish the conditions under which p c is a good approximation 
to the total perturbation pressure p. To this end consider the relation between p c , 
p 1 and p: 

^ 2 ||p C || 2 = lbH 2 + 7 2 MA||/|| 2 (47) 

which is justified if p(x) are p 7 (x) are assumed to be decorrelated. This is true by 
construction at t — 0 in our direct numerical simulation code (discussed in a later 
section). From Eq. (18), | |p 7 1 1 and ||u*|| are related by 

lb'll = ailb'ir, 

= «i(l - Xo), (48) 

where ai is an assumed 0(1) constant. Using Eq. (48), Eq. (47) becomes 

^ 2 |b c H 2 = lbo|| 2 + 7 2 MX(l-Xo) 2 - (49) 

The incompressible component of the pressure field is therefore negligible when 

I boll » 7^ai(l - Xo)- (50) 

This result will be verified with the help of direct numerical simulations. Equa- 
tion (50) tells us that one cannot neglect the influence of p l at higher levels of Mr 
and at lower levels of the initial compressibility ratio, Xo- 


III. Numerical Simulation 


A. Numerical Algorithm 

Our direct simulations of two-dimensional compressible turbulence are based on Eqs. (1) 
The flow is assumed to be isotropic. Hence a double Fourier representation is appropri- 
ate. Spectral methods, by nature of their high accuracy and low dispersive and dissi- 
pative errors are ideally suited for this problem. Spatial derivatives are approximated 
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by Fourier collocation (Canuto, Hussaini, Quarteroni and Zang 1988). In each coor- 
dinate direction, N grid points are used; for example Xj = 2-xj/N, j = 0, 1, • • • , N — 1 
in the x direction. The derivative of a function .F(x) with respect to x is calculated 
from the analytic derivative of the trigonometric interpolant of P(x) in the direction 
x. Many simulations of incompressible, homogeneous turbulence have used a Fourier- 
Galerkin method. The compressible equations, however, contain cubic rather than 
quadratic nonlinearities and true Galerkin methods are more expensive (compared 
with collocation methods) than they are for incompressible flow. The essential dif- 
ference between collocation and Galerkin methods is that the former are subject to 
both truncation and aliasing errors, whereas the latter have only truncation errors. 
As discussed extensively by Canuto, et al. (1988), the aliasing terms are not signifi- 
cant for a well-resolved flow. However, care is needed to pose the relevant equations 
for a collocation method in a form which ensures numerical stability. For this reason, 
the convective term in the momentum equations is actually used in the equivalent 
form 

-[V-(puu) + pu • Vu + uV-(pu)]. (51) 

As noted by Feiereisen, Reynolds and Ferziger (1981), this form, together with a 
symmetric differencing method in space (for example Fourier collocation), ensure 
conservation of mass and momentum. Moreover, energy is conserved for the semi- 
discrete inviscid equations. 


We are interested in the low subsonic regime Mr << 1 . The sound speed is in 
this case much greater than the flow velocity and this imposes a severe restriction 
on the time step of any explicit time marching numerical scheme. To remove this 
constraint, a splitting method is adopted. The first step integrates implicitly 


d(pu) 1 1 

—Qf- + rj" [V-(puu) + u • V (? u ) + P u • Vu] = — Va, (52) 

dP 

-qT + u • VP -(- 7 ?V u — CqV(pu) = 

— !— v 2 t + 7(7 ~ DMk $ 

RePr T Re 


from which the acoustic effects that vary on the fast time scale 0(M R ) have been 
removed. The second step integrates the “acoustic equations” which vary on the fast 
time scale. These equations are 


^ + -L_vp 


dt 

dP 

~di 


7 Ml 


+ CqV • (pu) 


0 , 

0 , 

0 . 


(53) 
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where c 0 is the root mean square (rms) of the sound speed. 


This splitting is employed at each stage of a low-storage third-order Runge-Kutta 
method. To minimize the errors due to the large implicit time step, Eqs. (53) are 
integrated analytically. In Fourier space the system (53) become 


dt 

dm 


+ tk • m 

ik 


+ P = 

dt 7 M R 

dP 2l „ 

-dt + ic ° k ’ m = 


0 , 

0 , 

0 , 


(54) 


where m = pu and Fourier transformed quantities, which depend upon the wavenum- 
ber k, are denoted by a circumflex. Equations (54) are solved exactly in Erlebacher, 
Hussaini, Speziale and Zang (1987) and are rewritten here for completeness: 


P W = 


i ^ r A ,CokAt s . £ , 


- ( 2 ) 

rrr ' 


p (2) 


— — 


zk * . c 0 kAt a * c 0 kAt, - 

l A sm ( ,, rr) ~ B COS ( 77 7= ) + B l 


c 0 kM R ^/y 


,cokAt a . * . ,C(jkAt a . 

A co <T1—7=) + B sm (Ti-^=)> 


'M Ry /i* 


(55) 




'M Ry / r 


where k = |k| is the magnitude of the Fourier wavevector, A = P , B — i^-k • m^ 1 ^. 
Superscripts 1 and 2 relate to the end of the first and second fractional step respec- 
tively. The effective time-step of the Runge-Kutta stage is denoted by A t a . The 
advantage of this splitting is that the principal terms responsible for the acoustic 
waves have been isolated. Since they are treated semi-implicitly, one expects the 
explicit time-step limitation to depend upon v + |c — c 0 | rather than v + c. (Although 
there is also a viscous stability limit for the first fractional step, it is well below the 
advection limit in the cases of interest.) This is clearly efficient at low Mach num- 
bers. Since the second fractional step is integrated analytically, it does not generate 
dispersion or dissipation errors. The only errors incurred are time splitting errors. 


If one is truly interested in the detailed time-evolution of the acoustic components, 
the time-step must be small enough to resolve the temporal evolution of these waves. 
Note however, that because the acoustics are treated exactly in the second step, the 
large time step simulation still produces accurate results, but the data is not available 
at intermediate times. 


The expected stability limit of this two-dimensional Fourier collocation method 
for the compressible Navier-Stokes equations has the form 


At < a min(- 


Ax 


+ 


Ay 


u| + |c - c 0 | M + |c - c 0 | 




(56) 
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For the third-order Runge-Kutta method employed here a = 0.54. For accuracy 
purposes, the results herein are all based ona = 0.1. 


B. Initial Conditions 

The initial conditions for the direct numerical simulations presented here are similar 
to those presented by Passot and Pouquet (1987). They are sufficiently general to 
produce the various turbulent regimes we expect to occur when Mr << 1. One must 
separately specify the reference Mach number Mr, the Reynolds number, Re, the 
Prandtl number, Pr, the rms levels of Uq , Uq, ||p 0 ||, ||T 0 ||, and the autocorrelation 
spectrum for p, T, u c and u 1 . We now detail the procedure used in the numerical 
code. 

We first generate random fields in physical space in the range [—0.5, 0.5] for u 0 , 
p 0 , and To. Next, the velocity is decomposed into the sum 

u o(x) = <(x) + u'(x) (57) 

where 


V-Uq = 0, 

Vxuq - 0. 


This decomposition is unique for homogeneous flows and is accomplished in Fourier 
space according to the prescription 


u 7 = 


u c = 


u 


k u 

k 2 


u — u 7 . 


After transforming the velocity components back to physical space, we rescale u c , 
u J , p, and T to impose a prescribed autocorrelation spectrum. The autocorrelation 
spectrum £pp(k) of a variable p(x) is defined as 

J p(x)p(x)dx = f Epp(k)dk. (58) 

For example, the 2D energy spectrum for u, i.e. E(k), is the autocorrelation spectrum 
E uu- 

We now impose the spectrum 

_ 

E(k) = k 4 e (59) 

on each of the variables. The procedure for imposing the required spectrum E(k) is 
explained using the density as an example. One scans the wavenumber space, and 
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for each pair (k x , k v ), the corresponding spherical shell is determined. The i th shell is 
defined by 

(t - i < k < i + i). (60) 

where k = + k^. All spherical shells have unit width. The 0 tfc shell is empty 

because k = (0,0) is the mean component which does not contribute to the autocor- 
relation spectrum. The contribution of the i th shell to the autocorrelation spectrum 
thus becomes 

e;= e wop (6i) 

where E* is the discrete representation of E*(k)dk. (For 2-D flows, spherical 

should be interpreted as circular.) Defining Ei as 


Ei = E(i) 

where E(i ) is obtained from Eq. (59), the density is rescaled according to 


(62) 


K') 



(63) 


in order to impose the desired spectrum. In a similar way, the autocorrelation spectra 
for u c , u 7 and T are calculated. For the velocity, the procedure is identical to that 
described for p, except that the scalar p is replaced by a vector, and |p(k)| 2 is replaced 
by|u(k)| 2 . 


It is still possible to rescale these quantities by an arbitrary constant in physical 
space without changing the properties of its autocorrelation spectrum. Therefore, we 
impose given rms levels for p, T and u by an appropriate rescaling. There still remains 
a degree of freedom on the velocity, since u c and u J can be weighted independently 
without affecting the rms of u. Therefore, we specify the initial level of compressibility 
in the flow according to 

_ / u c2 dx 

X / (u c2 + u /2 )dx 

The denominator of Eq. (64) is the total kinetic energy of the system since / u c -u 7 = 0 
by construction. 



The mean values of p and T are then added to the thermodynamic variables to 
obtain 


P < — 1 + Py 
T < — i + r. 


The reference pressure in the code is chosen to be Pr = prUr ; consequently the 
nondimensional pressure is then calculated from 

PT 

? 7MA 
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where p and T are respectively the non-dimensional density and temperature. The 
turbulent Mach number M t is related to the reference Mach number Mr according 
to Eq. (6). 


IV. Numerical Validation of Theory 


In this section, we show results of 2-D direct numerical simulations on grids of 64 x 
64 to validate the theory described in the preceding sections. This validation is 
accomplished by comparing the value of x after several acoustic periods as predicted 
by the computation, against that given by the theory. The parameter F, defined by 
Eq. (46) is also computed and compared with numerical results. 

Given the large parameter space, in every run, the density and temperature fluctu- 
ation rms levels are set equal to each other. An extensive parameter range is covered 
by considering all combinations of 

Xo = (0.1, 0.5, 0.9), 

Mr = (0.03,0.08,0.3), 

IWI-imi = (0.01,0.05,0.1). 

Cases are referred to by a 3 digit combination. For example, case 123 refers to 
Xo = 0.1, Mr = 0.08 and \\p\\ = ||T|| = 0.10. Simulations are run at Re = 150, 
||u|| = 1, Jfe 0 = 10. This corresponds to a microscale Reynolds number R\ — 20. The 
various diagnostics are computed at every iteration. 

The quantities considered for comparison between the numerical computations 
and the theory results are 



and x- Both these quantities are computed by averaging them over several consecutive 
iterations: x i fi averaged over iterations 5 to 10 (xs-io) and 10 to 15 (xio-is), while 
F is averaged over iterations 10-15 (Fio-is)- Their equilibrium values Xoo and Foo 
are computed by averaging their values over the last 60 iterations of the 400 iteration 
numerical simulation. The two sets of averages for x give an idea of the convergence 
rate of x towards the theoretical prediction x th • If X is averaged over too few time 
steps, the value could be either over or under predicted since x(0 is an oscillatory 
integral which tends to unity. On the other hand, if x is averaged over too many 
steps, the effects of viscosity will be noticeable, and agreement with x^ will not be 
reached. Averages over five iterations prove adequate. 

Table 3 summarizes the simulation results. Note that Xio -15 is in closer agreement 
with Xto u th an is Xb-io- On the other hand the average x decays as a function of 
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case 

X5-10 

XlO-15 

Xoo 

Aoo 

F 0 

-f 10-15 

F 

1 oo 

VIIPoll 

in 

0.11 

0.11 

0.08 

0.11 

1.00 

1.01 

1.01 

1.00 

112 

0.63 

0.62 

0.44 

0.62 

0.03 

0.98 

1.01 

1.00 

113 

0.90 

0.86 

0.76 

0.87 

0.01 

0.98 

0.99 

1.00 

121 

0.06 

0.06 

0.05 

0.06 

6.15 

1.00 

1.00 

1.05 

122 

0.25 

0.21 

0.14 

0.22 

0.25 

0.92 

1.00 

1.00 

123 

0.58 

0.52 

0.32 

0.49 

0.06 

1.68 

0.99 

1.00 

131 

0.08 

0.05 

0.06 

0.05 

86.5 

0.73 

0.93 

5.10 

132 

0.09 

0.07 

0.06 

0.07 

3.46 

0.88 

0.96 

1.39 

133 

0.11 

0.13 

0.09 

0.11 

0.87 

1.09 

0.95 

1.10 

211 

0.38 

0.38 

0.24 

0.38 

4.32 

1.03 

1.01 

1.00 

212 

0.77 

0.77 

0.59 

0.77 

0.17 

1.01 

1.00 

1.00 

213 

0.94 

0.92 

0.84 

0.92 

0.04 

0.96 

0.99 

1.00 

221 

0.30 

0.35 

0.21 

0.34 

30.7 

1.07 

1.01 

1.01 

222 

0.47 

0.47 

0.29 

0.47 

1.23 

0.93 

0.99 

1.00 

223 

0.73 

0.69 

0.29 

0.68 

0.31 

1.23 

1.00 

1 00 

231 

0.40 

0.27 

0.47 

0.33 

432. 

0.69 

0.96 

2.93 

232 

0.40 

0.29 

0.20 

0.35 

17.3 

0.71 

0.96 

1.12 

233 

0.42 

0.35 

0.23 

0.38 

4.33 

0.79 

0.99 

1.03 

311 

0.84 

0.84 

0.67 

0.83 

7.78 

1.03 

1.01 

1.00 

312 

0.95 

0.95 

0.88 

0.95 

0.31 

1.02 

1.00 

1.00 

313 

0.99 

0.98 

0.92 

0.98 

0.08 

0.93 

0.99 

1.00 

321 

0.78 

0.83 

0.63 

0.82 

55.3 

1.07 

1.03 

1.00 

322 

0.86 

0.86 

0.70 

0.87 

2.22 

0.96 

1.01 

1.00 

323 

0.94 

0.92 

0.81 

0.93 

0.55 

1.06 

0.99 

1.00 

331 

0.85 

0.76 

0.63 

0.82 

778. 

0.71 

1.02 

1.12 

332 

0.85 

0.77 

0.63 

0.82 

31.2 

0.72 

1.01 

1.00 

333 

0.86 

0.80 

0.64 

0.83 

7.80 

0.79 

1.01 

1.00 


Table 3: Comparison of two-dimensional direct simulation data against theo- 
retical predictions. All symbols are defined in the text. 
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time in part due to viscous effects which is reflected by the low values of x«> compared 
to Xtc table. 

Sarkar et al. (1989) showed analytically that, for low Mach number turbulence, 
the parameter F approaches, and then oscillates about a value of = 1. Table 
3, regarding F x are in agreement with this result. The value of F depends on the 
compressible component of pressure, 6p c , not p. Nonetheless, it is p, not p c which is 
used to compute F for table 3. Furthermore, the time interval, required to obtain a 
good average over the acoustic oscillations also increases with M^r. This explains why 
■^io-i5 is not close to unity at the higher values of Mr. The last column measures 
the departure of the initial compressible rms pressure from the initial rms pressure. 
The cases where £/||po|| > 1 correlate well with a non-equilibrium value of F after 15 
iterations. On the viscous time scale, F is nonetheless close to unity in all cases, which 
reflects that the turbulent Mach number has decreased substantially (see Eq. (50)), 
i.e. p 1 has become negligible with respect to p c . 

The ratio of Xio-is to x» is furthest from unity when M R = 0.3, which is explained 
by the longer acoustic time scale. Since the averages are computed over a fixed number 
of iterations, the sample length at Mr = 0.3 is not adequate, and the prediction for 
Xoo is not very good. The lower value of Xoo with respect to x£ arises because x 
decays on the viscous time-scale. On the other hand, F remains close to unity, and 
is basically unaffected by the viscosity terms. 

Turbulent simulations with initial flowfields for which F 0 = 1 preclude initial time 
boundary-layers which might otherwise impose stringent restrictions on time-stepping 
and code accuracy. We now derive an implicit relation for x as a function of Mr and 
IIpII- To this end, we begin by collecting the required formulas. Henceforth, all 
quantities refer to initial conditions that are consistent with an equilibrium ( F = 1). 
Based on a general velocity decomposition 

u = a(u c + pu 1 ) ( 67 ) 

where a is such that ||u|| = 1, the compressible to total kinetic energy ratio is written 
as 

x= ML_ = W 

I|u c II ! + /3 j I|u'|| ! 7 ! m£’ (68) 

and the rms of the compressible pressure field becomes 

lb c H 2 = IIpII 2 + 7 2 m*||p'|| 2 , 

= IIpII 2 

where we have assumed that p and p 1 are decorrelated. By construction, they are 
decorrelated at t — 0 in the numerical code from Eq. (21) More realistically how- 
ever, p 1 and p c are probably decorrelated since the incompressible and compressible 
components of the flowfield evolve quasi-independently. 
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From the previous set of equations one obtains 

M* = 7’A/£[x-a;M|(l-x) ! ]. ( 69 ) 

Contour plots of x are shown in Figure 1 (aq = 1). It is seen that for moderate to 
high values of x, the relationship between ||p|| and Mr is linear, which simply restates 
that ||p c || and ||p|| are approximately equal. At low values of x> ||p|| decreases, and 
eventually reaches zero. In other words, in a turbulent flow in acoustic equilibrium 
(F = 1), high rms pressure at high turbulent Mach numbers imply that % is c l° se 
to unity. Since a necessary condition for the generation of shocks is 8 — 0(1), it 
follows that when shocks are present in a turbulent compressible homogeneous flow 
in acoustic equilibrium, the compressible part of the flow must dominate. 


V. Conclusions 

In this paper, an asymptotic theory that is based on an initial- value formulation of the 
equations of motion for compressible homogeneous turbulent flow has been presented. 
The expansion parameter is the turbulent Mach number, which is assumed to be 
small. A key feature of the theory is the decomposition of the turbulent velocity field 
into two physically meaningful components. The first is solenoidal and satisfies the 
time-dependent incompressible Navier-Stokes equations. This uniquely determines 
the second component which is initially irrotational, but acquires a small amount of 
vorticity at later times. One consequence of this decomposition, is the simultaneous 
separation of the pressure field into incompressible and compressible components. 
The incompressible component is proportional to M 

To lowest order, the asymptotic theory describes the acoustic part of the flow by a 
wave equation which is coupled to the solenoidal flow through the initial conditions. 
At later times, stronger coupling will occur through the neglected convective and 
viscous terms in the momentum equation. Order of magnitude estimates permit an 
apriori estimate of the compressibility levels of the flow as a function of the initial 
conditions. It is found that the initial flow can evolve towards two basically different 
states (in the absence of shocks). If Fq < 1, the pressure fluctuations remain at their 
initial levels, while the velocity fluctuations increase in such a way that x becomes 
0(1). On the other hand, if F 0 > 1, large pressure waves develop over a time of 
0(Mr)j and the rms level of the acoustic velocity remains constant. This transient 
behavior is well described by the solution to the wave equation, and has been checked 
against extensive two-dimensional direct numerical simulations. Note however that 
the theory is also valid for three-dimensional turbulent flows and spot checks have 
been conducted. 

It has also been established that a necessary, but not sufficient, condition for the 
generation of shocks is that the flow must be in an initial state of disequilibrium (i.e. 
F < 1) which results from an initial pressure level which is 0(1). 
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Several issues need further clarification. First, it is still not clear if there are more 
stringent initial conditions which could guarantee the emergence of isotropic shock 
distributions in two or three-dimensional compressible turbulent flows. Moreover, 
it is not clear how the theory should be modified if the reference Mach number is 
no longer small. Second, the current paper only treated the lowest order equations 
which resulted from the asymptotic expansions. To understand the evolution of the 
turbulent flow on the convective and viscous time scales, it is necessary to keep both 
viscous terms, and higher order terms in the Mach number expansion. These higher 
order terms will allow an explicit computation of the effects of the incompressible 
eddies on the acoustic waves and vice versa. These two issues, are now being studied 
and will be the focus of future publications. 
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Appendix A 


I. Appendix 


We want to explain the behavior of hyperbolic systems which are far from symmetric. 
Let us first consider the system of ordinary differential equations 

ui= ' 4u= (- ( l o )"■ u= (“0’ <ai) 

with initial data 

u(0) = u 0 . (A2) 

The matrix A is antisymmetric ( A T = —A) so the eigenvalues fij, (j = 1, 2) are purely 
imaginary and the corresponding eigenvectors y } are orthonormal with respect to the 
scalar product 


< u, v >= UjU! + u 2 v 2 , |u| 2 =< u, u >, (A3) 

where an overbar denotes complex conjugate. The general solution of Eq. (Al) is 
given by 

u(t) = a 1 e tfilt y l + a 2 e Lfllt y 2 (A4) 

where 

o-i=<yi»uo>, cr 2 =< y 2 , u 0 > . (A5) 

From the orthogonality of the eigenvectors, 


l u (0l 2 — kil 2 + l°2| 2 = l u (0)| 2 > 


(A6) 


and the amplitude of the solution does not change. This is true for any antisymmetric 
matrix which can be shown by the energy method 



d 

a? < u,u > 

< u e ,u > + < u,u t > 

< Au,u > + < u,4u > 

— < u,Au > + < u,Au > 

0 . 
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By the previous result, 

|u(t )| 2 = |ui(t )| 2 + \u2(t )\ 2 

= M 0)| 2 + |n 2 ( 0)| 2 

= |m(0)| a + e|ii a (0)| 2 . 

Assume now that u^O) = 0(1) and u 2 (0) = 0(1). The general form (A4) of the 
solution tells us that in general both Uj(t) = 0(1) and ^(t) — 0(1). Therefore we 
obtain for the original variables 

ui(t) = ui(t) - 0(1) 
u 2 (t) = e~^ 2 u 2 (t) = oor 1 / 2 ) 

which shows that the amplitude 

KOI = °( e " 1 /! )I u (o)I (aio) 

increases dramatically for sufficiently small e though the eigenvalues of A are still 
purely imaginary. We can also express the result in the following way. Assume that 
u x ( 0 ) = 0 and u a (0) = 1. Then |u(i)| 2 = e and 

KOI = (All) 

i.e. for this initial data, the amplitude does not increase. If we now make a small 
perturbation in ui( 0 ), say ui( 0 ) = 77 , then the perturbation will be amplified such that 
m(t ) ~ 0(1 + e -1 / 2 ?;) and v^(t) ss 0(1 + e - 1 / 2 7 7 ). Thus if 77 >> e 1 / 2 the perturbation 
destroys the solution. 

There are no difficulties in generalizing the results. Consider the system 

u t + Au = 0, (A12) 
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where A is an n x n matrix and u is a vector with n components. Assume that there 
is a diagonal transformation 

D — diag(di, • • • , d n ), 0 < d, < 1, (M3) 

such that DAD -1 is antisymmetric. Then the amplitude of the solution can grow by 
a factor maXj d“ , i.e. 

| u (t)| = 0(m|.xdj 1 ) |u(0)|. (A14) 

In particular, if we initialize the data such that 

|u(f)| « |u(0)|, (A15) 

then we can find perturbations rj which grow to the order maXj dj 1 ?]. 

After Fourier transform, a hyperbolic system 



The Fourier transform is of the form (A7). 
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Figure 1: ||p(M^, %)|| as a function of Mr for different values of x* The contour levels 
of x all correspond to the condition that F = 1. 
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